#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*
 * 构造拟合多项式
 * 参数一：散点数组，第二维下标0指X值，下标1指Y值
 * 参数二：散点个数
 * 参数三：所指定的拟合次数
 * 参数四：为传出参数，传出构造的拟合多项式的平方误差
 * 返回值：构造的对应的拟合多项式数组首地址，下标即为对应项次数
 */
double *ConstructFittingPolynomials(double Scatter[][2],int NumOFScatter,int NumberOfFitting,double *SquareError);

//普通高斯消元法通解
double* GaussianEliminationGEN(double **Matrix,int line,int column);
//打印矩阵
void printMatrix(double **Matrix,int line,int column);

//打印多项式
void printPolynomial(double Polynomial[],int N);

int main() {

    //测试样例一：通过
//    double Scatter1[7][2] = {{19.1,76.3},{25,77.8},{30.1,79.25},{36,80.8},{40,82.35},{45.1,83.9},{50,85.1}};
//    double *Ans;
//    Ans = ConstructFittingPolynomials(Scatter1,7,1);
//    printPolynomial(Ans,2);

    //测试样例二：通过
//    double Scatter2[9][2] = {{1,10},{3,5},{4,4},{5,2},{6,1},{7,1},{8,2},{9,3},{10,4}};
//    double *Ans;
//    Ans = ConstructFittingPolynomials(Scatter2,9,2);
//    printPolynomial(Ans,3);

    //例题：
//    double Scatter2[7][2] = {{0,1},{0.5,1.75},{0.6,1.96},{0.7,2.19},{0.8,2.44},{0.9,2.71},{1,3}};
//    double *Ans,SquareError;
//    Ans = ConstructFittingPolynomials(Scatter2,7,2,&SquareError);
//    printPolynomial(Ans,3);
//    printf("Square Error is :%lf\n",SquareError);

    int ScatterNUM,FittingNUM;  //ScatterNUM是输入的散点个数，FittingNUM是要进行的拟合次数
    printf("please input Scatter NUM And Fitting NUM:\n");
    scanf("%d%d",&ScatterNUM,&FittingNUM);
    double Scatter[ScatterNUM][2];      //用于存放输入的散点
    printf("please input Scatter point: X AND Y\n");
    for (int i = 0; i < ScatterNUM; ++i)
        scanf("%lf%lf",&Scatter[i][0],&Scatter[i][1]);
    double *Ans,SquareERR;  //Ans用于接受构造的拟合多项式数组，SquareERR用于存放拟合多项式的平方误差
    Ans = ConstructFittingPolynomials(Scatter,ScatterNUM,FittingNUM,&SquareERR);
    printPolynomial(Ans,FittingNUM+1);  //将拟合多项式打印出来
    printf("Square Error is :%lf\n",SquareERR); //输出平方误差

    free(Ans);
    return 0;
}

double *ConstructFittingPolynomials(double Scatter[][2],int NumOFScatter,int NumberOfFitting,double *SquareError){
    int XN_NUM = 2*(NumberOfFitting+1)-1,YN_NUM = NumberOfFitting+1;
    double XN[XN_NUM];   //用于存放左边矩阵中各个X^n求和的值
    double YN[YN_NUM];   //用于存放右边矩阵中各个X*Y求和的值
    double *FittingPolynomial ;    //用于存储最终拟合多项式
    memset(XN,0, sizeof(double )*2*NumberOfFitting);    //初始化值
    memset(YN,0, sizeof(double )*(NumberOfFitting+1));

    //填补XN，YN
    for (int i = 0; i < XN_NUM; ++i) {
        for (int j = 0; j < NumOFScatter; ++j) {
            XN[i] += pow(Scatter[j][0],i);
            if(i < YN_NUM)
                YN[i] += Scatter[j][1]* pow(Scatter[j][0],i);
        }
    }

    //构造增广矩阵
    //手动构建二维增广矩阵
    double **FittingAugmentedMatrix = (double **)malloc(sizeof(double *)*(NumberOfFitting+1));
    for (int i = 0; i < NumberOfFitting+1; ++i) {
        FittingAugmentedMatrix[i] = (double *) malloc(sizeof(double )*(NumberOfFitting+2));
    }
    //填补增广矩阵
    for (int i = 0; i < NumberOfFitting+1; ++i) {
        for (int j = 0,k=i; j < NumberOfFitting+1; ++j,k++)
            FittingAugmentedMatrix[i][j] = XN[k];
        FittingAugmentedMatrix[i][NumberOfFitting+1] = YN[i];
    }

    //高斯消元法解增广矩阵求拟合多项式的各项系数
    FittingPolynomial = GaussianEliminationGEN(FittingAugmentedMatrix,NumberOfFitting+1,NumberOfFitting+2);

    //求平方误差
    *SquareError = 0;
    for (int i = 0; i < NumOFScatter; ++i) {
        double px=0;
        for (int j = 0; j < NumberOfFitting+1; ++j) {
            px += FittingPolynomial[j]* pow(Scatter[i][0],j);
        }
        *SquareError += fabs(px-Scatter[i][1])* fabs(px-Scatter[i][1]);
    }

    //释放堆内存
    for (int i = 0; i < NumberOfFitting+1; ++i)
        free(FittingAugmentedMatrix[i]);
    free(FittingAugmentedMatrix);

    return FittingPolynomial;   //返回多项式
}

void printPolynomial(double Polynomial[],int N){
    for (int i = 0; i < N; ++i) {
        printf("(%.5lf)x^%d\t",Polynomial[i],i);
    }
    printf("\n");
}

//打印矩阵
void printMatrix(double **Matrix,int line,int column){
    for (int i = 0; i < line; ++i) {
        for (int j = 0; j < column; ++j) {
            printf("%lf\t",Matrix[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

//普通高斯消元法通解
double* GaussianEliminationGEN(double **Matrix,int line,int column){
    double key,*Ans;  //Ans是用于存放最终答案的数组首地址，key是高斯消元消每列当前行的因子，即当前列当前行除以当前列的主元素的值
    Ans = malloc(sizeof(double )*line); //用于给答案存放开辟空间

    for (int i = 1; i <= column-2; ++i) {   //一共需要消除的列数
        for (int j = i+1; j <= line; ++j) {     //每一列需要消除的行数，j为目前所需要消除的行号

            printMatrix(Matrix,line,column);    //消除之前打印一下当前矩阵状况

            key=-Matrix[j-1][i-1]/Matrix[i-1][i-1]; //获取消除每行当前列的消除因子
            for (int k = 0; k < column; ++k)    //对当前行做变换
                Matrix[j-1][k]+=key*Matrix[i-1][k];
        }
    }
    printMatrix(Matrix,line,column);    //变换完成之后再次打印矩阵，查看当前矩阵状况

    for (int i = 0; i < line; ++i) {    //总共需要计算的答案数
        double sum=Matrix[line-1-i][column-1];  //回代答案被当前求的未知数的系数除之前总和的初值
        for (int j = line-1; j > line-1-i ; --j) {  //对每行答案回代并且右移，最终获取总和的终值
            sum-=Matrix[line-1-i][j]*Ans[j];
        }
        Ans[line-1-i]=sum/Matrix[line-1-i][line-1-i];   //获取当前所求未知数的值，即总和除以所求未知数系数
    }

    return Ans; //返回答案
}